
% A2c_price_index_full_adj
%==========================================================================

% Description: This code computes product-specific price indexes in the
% full model, both factually and counterfactually 

clear
clc

cd 'E:\data_replication'

% True Utilities
%==========================================================================

% Import Demand Parameter Estimates
%--------------------------------------------------------------------------

load('estimation/5_supply_side/main_data_supply.mat');
clearvars -except i_k index_price_coeff X_MPEC_opt_summary Y_summary K K_A

for i_k = 1:4092
   
i_k 

filename_costs = sprintf('robustness/firm_heterogeneity/cost_estimates/costs_%i_k.mat', i_k);  
filename_costs_2 = sprintf('robustness/firm_heterogeneity/cost_estimates/costs_full_adj_%i_k.mat', i_k);  

if exist(filename_costs, 'file') == 2
if exist(filename_costs_2, 'file') == 2
    
k = index_price_coeff(i_k,1);    
X_MPEC_opt = X_MPEC_opt_summary{k};
Y_k = Y_summary{k};
alpha_constant = X_MPEC_opt(1,1);
alpha_distance = X_MPEC_opt(2,1);
alpha_population = X_MPEC_opt(3,1);
alpha_price = X_MPEC_opt(4,1);
alpha_1 = X_MPEC_opt(12,1);
alpha_FE = X_MPEC_opt((3*(K_A+1) + 1):(3*(K_A+1) + (K - K_A)), 1);


% Import Cost Estimates
%--------------------------------------------------------------------------

filename_costs_matlab = sprintf('robustness/firm_heterogeneity/cost_estimates/costs_%i_k.mat', i_k);  
load(filename_costs_matlab)

numProdsTotal = size(data_save,1);

data_k = [data_save.Year, data_save.Quarter, data_save.Declarant];
declarant = data_save.Declarant;
price_k = data_save.price_k;
q_k = data_save.q_k;
share_k = data_save.share_k;
distance_k = data_save.distance_k;
population_log_k = data_save.population_log_k;
mc_flag_true = data_save.p_flag;
Y_k = Y_summary{k};
market_id_k = data_save.market_id;
T = market_id_k(end,1);
nn = size(Y_k,2);


prodsMarket = zeros(T,1);
for t = 1:T
prodsMarket_temp = sum(market_id_k==t);
prodsMarket(t,1) = prodsMarket_temp(1,1);
end

marketStarts = zeros(T,1);
marketEnds = zeros(T,1);
marketStarts(1) = 1;
marketEnds(1) = prodsMarket(1);
for t=2:T
    marketStarts(t) = marketEnds(t-1) + 1;
    marketEnds(t) = marketStarts(t) + prodsMarket(t) - 1;
end

marketForProducts = zeros(numProdsTotal,1);
for t=1:T
marketForProducts(marketStarts(t):marketEnds(t)) = t;
end

for t=1:T
    data_market = [data_save.Year, data_save.Quarter, data_save.Declarant, data_save.Partner];
    data_market = data_market(marketStarts(t):marketEnds(t), :);
    N_true(t,1) = size(data_market, 1);    
    N_dom_true(t,1) = sum(data_market(:,3) == data_market(:,4));
end

delta_k =  q_k + log(price_k)*alpha_price + alpha_population*population_log_k + alpha_distance*distance_k;

for t = 1:T   
Y_market = Y_k(t,:);
expmu(marketStarts(t):marketEnds(t),:) = exp(log(price_k(marketStarts(t):marketEnds(t),:))*alpha_1*Y_market);
end
expmeanval = exp(delta_k);
oo = ones(1,nn);                  
sharesum = sparse(zeros(T,numProdsTotal));  % used to create denominators in logit predicted shares (i.e. sums numerators)
for t = 1:T
    sharesum(t,marketStarts(t):marketEnds(t)) = 1;
end
numer = (expmeanval*oo ).*expmu;        
sum1 = sharesum*numer;
sum11 = 1./(0+sum1);                                                       % Will not result in share_k = EstShare_true
denom1 = sum11(marketForProducts,:);    
simShare_true = numer.*denom1;              
EstShare_true = mean(simShare_true,2);  

v_true = log(numer); 
P_true = zeros(T,nn);
beta_price = Y_k * alpha_1 + alpha_price;
for tt = 1:T
v_true_t = v_true(marketStarts(tt,1):marketEnds(tt,1),:,:); 
beta_price_t = beta_price(tt, :);
exp_v_true_t = exp(v_true_t);
s_exp_v_true_t = sum(exp_v_true_t, 1);
ln_s_exp_v_true_t = log(s_exp_v_true_t);
P_true(tt,:) = exp(1) * exp(1./beta_price_t .* ln_s_exp_v_true_t) .* gamma(1 - 1./beta_price_t);
end

mc_flag_true = mc_flag_true(marketStarts, 1);

clearvars -except i_k P_true mc_flag_true N_true N_dom_true 


% Import f_opt flag
%==========================================================================

filename_f_opt_matlab = sprintf('robustness/firm_heterogeneity/cost_estimates/data_f_opt_%i_k.mat', i_k);  
load(filename_f_opt_matlab);
mc_flag_f_opt = data_f_opt.p_flag;


% Counterfactual Utilities
%==========================================================================

load('estimation/5_supply_side/main_data_supply.mat');
clearvars -except i_k index_price_coeff X_MPEC_opt_summary Y_summary K K_A P_true mc_flag_true mc_flag_f_opt N_true N_dom_true 

k = index_price_coeff(i_k,1);    
X_MPEC_opt = X_MPEC_opt_summary{k};
Y_k = Y_summary{k};
alpha_constant = X_MPEC_opt(1,1);
alpha_distance = X_MPEC_opt(2,1);
alpha_population = X_MPEC_opt(3,1);
alpha_price = X_MPEC_opt(4,1);
alpha_1 = X_MPEC_opt(12,1);
alpha_FE = X_MPEC_opt((3*(K_A+1) + 1):(3*(K_A+1) + (K - K_A)), 1);


filename_costs_matlab = sprintf('robustness/firm_heterogeneity/cost_estimates/costs_full_adj_%i_k.mat', i_k);  
load(filename_costs_matlab)

numProdsTotal = size(data_save,1);

data_k = [data_save.Year, data_save.Quarter, data_save.Declarant];
declarant = data_save.Declarant;
partner = data_save.Partner;
price_k = data_save.price_k;
q_k = data_save.q_k;
share_k = data_save.share_k;
distance_k = data_save.distance_k;
population_log_k = data_save.population_log_k;
mc_flag_counter = data_save.p_flag;
Y_k = Y_summary{k};
market_id_k = data_save.market_id;
T = market_id_k(end,1);
nn = size(Y_k,2);


prodsMarket = zeros(T,1);
for t = 1:T
prodsMarket_temp = sum(market_id_k==t);
prodsMarket(t,1) = prodsMarket_temp(1,1);
end

marketStarts = zeros(T,1);
marketEnds = zeros(T,1);
marketStarts(1) = 1;
marketEnds(1) = prodsMarket(1);
for t=2:T
    marketStarts(t) = marketEnds(t-1) + 1;
    marketEnds(t) = marketStarts(t) + prodsMarket(t) - 1;
end

marketForProducts = zeros(numProdsTotal,1);
for t=1:T
marketForProducts(marketStarts(t):marketEnds(t)) = t;
end

for t=1:T
    data_market = [data_save.Year, data_save.Quarter, data_save.Declarant, data_save.Partner];
    data_market = data_market(marketStarts(t):marketEnds(t), :);
    N_counter(t,1) = size(data_market, 1);    
    N_dom_counter(t,1) = sum(data_market(:,3) == data_market(:,4));
end

delta_k =  q_k + log(price_k)*alpha_price + alpha_population*population_log_k + alpha_distance*distance_k;

for t = 1:T   
Y_market = Y_k(t,:);
expmu(marketStarts(t):marketEnds(t),:) = exp(log(price_k(marketStarts(t):marketEnds(t),:))*alpha_1*Y_market);
end
expmeanval = exp(delta_k);
oo = ones(1,nn);                  
sharesum = sparse(zeros(T,numProdsTotal));  % used to create denominators in logit predicted shares (i.e. sums numerators)
for t = 1:T
    sharesum(t,marketStarts(t):marketEnds(t)) = 1;
end
numer = (expmeanval*oo ).*expmu;        
sum1 = sharesum*numer;
sum11 = 1./(0+sum1);                                                       % Will not result in share_k = EstShare_true
denom1 = sum11(marketForProducts,:);    
simShare_true = numer.*denom1;              
EstShare_true = mean(simShare_true,2);  


v_counter = log(numer); 
P_counter = zeros(T,nn);
beta_price = Y_k * alpha_1 + alpha_price;
for tt = 1:T
v_counter_t = v_counter(marketStarts(tt,1):marketEnds(tt,1),:,:); 
beta_price_t = beta_price(tt, :);
exp_v_counter_t = exp(v_counter_t);
s_exp_v_counter_t = sum(exp_v_counter_t, 1);
ln_s_exp_v_counter_t = log(s_exp_v_counter_t);
P_counter(tt,:) = exp(1) * exp(1./beta_price_t .* ln_s_exp_v_counter_t) .* gamma(1 - 1./beta_price_t);
end

mc_flag_counter = mc_flag_counter(marketStarts, 1);

% Sort results
%--------------------------------------------------------------------------

beta_price = Y_k * alpha_1 + alpha_price;
ratio = P_counter ./ P_true;

ratio_sorted = zeros(T, nn);
P_true_sorted = zeros(T, nn);
P_counter_sorted = zeros(T, nn);
beta_price_sorted = zeros(T, nn);
for tt = 1:T
[Y_k_tt, ind] = sort(Y_k(tt,:));
P_true_sorted(tt,:) = P_true(tt,ind);
P_counter_sorted(tt,:) = P_counter(tt,ind);
ratio_sorted(tt,:) = ratio(tt,ind);
beta_price_sorted(tt,:) = beta_price(tt,ind);
end


diff = ratio_sorted(:,20) - ratio_sorted(:,80);
P_true_sorted_10 = P_true_sorted(:,10);
P_true_sorted_20 = P_true_sorted(:,20);
P_true_sorted_80 = P_true_sorted(:,80);
P_true_sorted_90 = P_true_sorted(:,90);
P_counter_sorted_10 = P_counter_sorted(:,10);
P_counter_sorted_20 = P_counter_sorted(:,20);
P_counter_sorted_80 = P_counter_sorted(:,80);
P_counter_sorted_90 = P_counter_sorted(:,90);
beta_price_20 = beta_price_sorted(:,20);
beta_price_80 = beta_price_sorted(:,80);

%==========================================================================


% Save Output to file
%==========================================================================

data_summary = zeros(T,3);
domestic_producer = ones(T,1);
for tt = 1:T
data_tt = data_k(marketStarts(tt,1):marketEnds(tt,1),:);     
data_summary(tt,:) = data_tt(1,:);
partner_tt = partner(marketStarts(tt,1):marketEnds(tt,1),1);  
if partner_tt(1,1) == 9999
domestic_producer(tt,1) = 0;
end
end


P_true_sorted_20_real = P_true_sorted_20 .* (imag(P_true_sorted_20) == 0);
P_true_sorted_20_real(isnan(P_true_sorted_20)) = 0;
P_true_sorted_20_real(P_true_sorted_20_real == Inf) = 0;
P_true_sorted_20_real(P_true_sorted_20_real == -Inf) = 0;
P_true_sorted_20_real(isnan(P_true_sorted_20_real)) = 0;
P_true_sorted_80_real = P_true_sorted_80 .* (imag(P_true_sorted_80) == 0);
P_true_sorted_80_real(isnan(P_true_sorted_80)) = 0;
P_true_sorted_80_real(P_true_sorted_80_real == Inf) = 0;
P_true_sorted_80_real(P_true_sorted_80_real == -Inf) = 0;
P_true_sorted_80_real(isnan(P_true_sorted_80_real)) = 0;

P_counter_sorted_20_real = P_counter_sorted_20 .* (imag(P_counter_sorted_20) == 0);
P_counter_sorted_20_real(isnan(P_counter_sorted_20)) = 0;
P_counter_sorted_20_real(P_counter_sorted_20_real == Inf) = 0;
P_counter_sorted_20_real(P_counter_sorted_20_real == -Inf) = 0;
P_counter_sorted_20_real(isnan(P_counter_sorted_20_real)) = 0;
P_counter_sorted_80_real = P_counter_sorted_80 .* (imag(P_counter_sorted_80) == 0);
P_counter_sorted_80_real(isnan(P_counter_sorted_80)) = 0;
P_counter_sorted_80_real(P_counter_sorted_80_real == Inf) = 0;
P_counter_sorted_80_real(P_counter_sorted_80_real == -Inf) = 0;
P_counter_sorted_80_real(isnan(P_counter_sorted_80_real)) = 0;

P_true_sorted_20 = P_true_sorted_20_real;
P_true_sorted_80 = P_true_sorted_80_real;
P_counter_sorted_20 = P_counter_sorted_20_real;
P_counter_sorted_80 = P_counter_sorted_80_real;


filename_table = sprintf('robustness/firm_heterogeneity/price_index/product_specific/diff_full_adj_table_%i_k.csv', i_k); 
data_save = array2table([data_summary, P_true_sorted_10, P_true_sorted_20, P_true_sorted_80, P_true_sorted_90, P_counter_sorted_10, P_counter_sorted_20, P_counter_sorted_80, P_counter_sorted_90, beta_price_20, beta_price_80, diff, domestic_producer, mc_flag_true, mc_flag_counter, mc_flag_f_opt, N_true, N_dom_true, N_counter, N_dom_counter]);
data_save.Properties.VariableNames = {'Year', 'Quarter', 'Declarant', 'P_true_sorted_10', 'P_true_sorted_20', 'P_true_sorted_80', 'P_true_sorted_90', 'P_counter_sorted_10', 'P_counter_sorted_20', 'P_counter_sorted_80', 'P_counter_sorted_90', 'beta_price_20', 'beta_price_80', 'diff', 'domestic_producer', 'mc_flag_true', 'mc_flag_counter', 'mc_flag_f_opt', 'N_true', 'N_dom_true', 'N_counter', 'N_dom_counter'};  
writetable(data_save, filename_table, 'Delimiter',',','QuoteStrings',true);

clearvars -except Y_summary index_price_coeff data_summary FE_summary i_k income_mu_sigma  K K_A KK objective_MPEC_summary price_coeff_MPEC price_table status_MPEC_index status_MPEC_summary v_summary X_MPEC_opt_summary X_MPEC_summary

end
end
end

